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We implement a Wang-Landau sampling technique in quantum Monte Carlo (QMC) for the pur- 
pose of calculating the Renyi entanglement entropies and associated mutual information. The algo- 
rithm converges an estimate for an analogue to the density of states for Stochastic Series Expansion 
QMC allowing a direct calculation of Renyi entropies without explicit thermodynamic integration. 
We benchmark results for the mutual information on two-dimensional (2D) isotropic and anisotropic 
Heisenberg models, 2D transverse field Ising model, and 3D Heisenberg model, confirming a critical 
scaling of the mutual information in cases with a finite-temperature transition. We discuss the 
benefits and limitations of broad sampling techniques compared to standard importance sampling 
methods. 



I. INTRODUCTION 

Entanglement is a measure of the ties that bind quan- 
tum mechanical particles across space and time. In a 
quantum many-body wavefunction, the amount of en- 
tanglement between two sub-regions of a system, A and 
B, is often quantified through a functional weighting of 
the eigenvalues of the reduced density matrix pA ~ 
TtbP, where there is always some set of such that 
p = K\4'i){4>i\- This is used to define the so-called 
Renyi entanglement entropies [1], 
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where the limit a = 1 is the familiar von Neumann en- 
tanglement entropy. The scaling of 5*1 in particular has 
been used to study and classify the groundstate phases 
and phase transitions of condensed matter systems 
[7]. Until recently it could only be accessed in numerical 
simulations through complete knowledge of the ground- 
state wavefunction - restricting its study to Hamiltonians 
which can be solved by exact diagonalization or related 
methods (e.g. density matrix renormalization group) [8]. 

In 2010, Hastings et. al. introduced a procedure to 
calculate Renyi entropies for integer a > 2 in zero- 
temperature projector quantum Monte Carlo (QMC) 
simulations through the expectation value of a swap op- 
erator [3]. This operator acts on two independent copies 
of the QMC-sampled configuration, literally "swapping" 
basis states in region A between the two copies. Since 
introduction of the algorithm, there has been a flurry 
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of activity in calculating Renyi entropies numerically us- 
ing swap- related techniques [TUlfH] . Soon thereafter [13] , 
it also was demonstrated that Renyi entropies could be 
calculated in finite-temperature QMC through adapta- 
tion of a well known "replica" trick [16l [17] , whereby the 
d+1 dimensional QMC simulation cell is doubled (in the 
case of S2) in size in its imaginary time direction [181119) . 
The Renyi entropy is then related to the logarithm of the 
ratio of the partition function of the replicated system 
{Z[A,2,T]) to the partition function of the original sys- 
tem, Z{T). Unlike the swap method that gives a Renyi 
entropy directly at T = 0, the replica trick calculated for 
two simulations (one for Z[A,2,Tq] and one for Z{Tq)) 
alone does not give the Renyi entropy at the tempera- 
ture Tq. Typically, a careful integration over many sim- 
ulations from T = 00 to T = Tq must be performed to 
calculate the Renyi entropy [TH [50] . Other approaches 
are possible, including a recent idea that allows direct 
calculation of the ratio of two partition functions at a 
fixed temperature without integration [21j . which effec- 
tively implements ideas of the swap operator at finite 
temperature. Ultimately, such ideas should be wed into 
the application of broad histogram techniques, which are 
natural solutions for Monte Carlo simulations which re- 
quire calculating the partition function. 

In this paper, we use an advanced Wang-Landau [^ 
sampling technique for calculating the Renyi entropy in 
finite-temperature QMC that avoids the need for many 
separate simulations, and allows access to the Renyi en- 
tropy at any arbitrary temperature between T = 00 and 
some chosen cut-off temperature. We outline the tech- 
nique for the specific case of S2 using a Stochastic Se- 
ries Expansion (SSE) QMC. We then demonstrate the 
method by looking at the finite-size scaling of the Renyi 
mutual information in three variants of the spin- 1/2 XXZ 
model in two and three dimensions and the 2D transverse 
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field Ising model. 



II. METHODS 



In this section we discuss the specific implementation 
details of the Wang-Landau method in the context of 
QMC and calculating the mutual information between 
two sub-regions A and B. This particular implementa- 
tion is guided by the estabhshed SSE QMC p3ti25] and 
recent extensions using a modified simulation topology 
[T5] . For adding Wang-Landau sampling to the SSE 
QMC we find the method outlined by Troyer, Wessel 
and Alet (2003) [26j is a simple and concise starting 
point, although it has similar limitations that have been 
found also to exist in classical Wang-Landau [57]. Go- 
ing beyond Wang-Landau approaches to more general 
broad histogram sampling, we implement a more ad- 
vanced technique detailed by Wessel, Stoop, Gull, Trebst, 
and Troyer (2007) [28^ that addresses some of the short- 
comings with the original Wang-Landau approach, but 
still provides the density of states needed to calculate 
thermodynamic properties (including the Renyi entropy) 
over a large range of temperatures. 



A. Stochastic series expansion 

The SSE [211 m] is a framework similar conceptually 
to world line methods [30H32] in which one simulates a d- 
dimensional quantum system using a {d + l)-dimensional 
classical simulation. To begin, we recall the standard SSE 
implementation, which represents the partition function 
as 



Z = Tr 
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where represents the trace over some complete basis 
of states for the system. We rewrite the above by collect- 
ing all the prefactors to the left and inserting the reso- 
lution of the identity between each power of the Hamil- 
tonian. Since the Hamiltonian is a sum, we can group 
operators into pieces that act on certain groups of sites 
and by whether the are diagonal or off diagonal in our 
basis of choice. These grouped pieces of the Hamiltonian 
we refer to as "bond" operators Ha^b where a differen- 
tiates diagonal and off diagonal terms and h denotes a 
grouping of sites. For our formalism, a = 1 represents 
diagonal terms while a — 2 represents off diagonal terms, 
while b specifies the bond connecting two sites that the 
Hamiltonian piece operates on. Using this formalism we 
obtain the SSE representation of the partition function. 
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where I'epresents the sum over all valid configura- 

tions {■00 • • ■ V'ra+i} with the condition (^Voj = (V'n+i . If 
we allow the sum to go from zero to infinity and include 
all possible fists of operators Ha^^bi for each n, this for- 
mulation is exact. In practice, contributions above some 
large value N are numerically insignificant and can be ig- 
nored. The cutoff N is always larger than any visited n 
in the importance sampling SSE QMC, and as such does 
not affect the statistics except to allow us to represent 
the set of operators as a finite list. 

The two updates we use to efficiently sample only non- 
zero weights are referred to as the diagonal update and the 
directed loop update j 24j . In the diagonal update we insert 
and remove diagonal operators from the list, changing n 
and the total weight, but not changing (V'o|- We use the 
change in weight at each insertion or removal and the 
usual Boltzmann condition to accept or reject these pro- 
posed moves. In the directed loop updates we connect 
all the Hamiltonian elements by what sites they interact 
with at each insertion of the identity in the expanded 
trace We then allow a loop to traverse this linked 

network of sites and change the operators (and spin con- 
figurations) it passes through, changing the type of each 
operator (and hence total weight) but not changing which 
sites each operator interacts with. Using these two up- 
dates our simulation is able to satisfy both ergodicity and 
detailed balance. 

Recently, the SSE QMC was adapted to allow for the 
calculation of Renyi entropy between two regions A and 
its complement B [20j. To do this, the so-called "replica 
trick" [33] was employed to estimate Eq. ([T]): 

Z{A,a,l3) 
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allowing the calculation the Renyi entropies through: 
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{pa) =. (ln[Z(A, a, /?)] - a ln[Z(/3)]) , (5) 

Here, the two partition functions represent two separate 
QMC simulations: Z{P) is the "usual" partition function 
for the physical system under study, and Z{A^ a, /?) is a 
replicated partition function of the form, 



Z(A,«,T) = >:(0o^|(^„^|e-''^|0f)|V'„^) 



(6) 



where ("0/^1 represents the basis containing spins in A, 
while I represent the spins not in region A and, as il- 
lustrated in Fig.[l] the boundary of region A is "stitched" 
together between the two replicas in the imaginary time 
direction. 

It is apparent from Eq. ^ that the calculation of Renyi 
entropies involves a procedure very similar to calculating 
a difference of free energies. As is well know, neither 
the free energy, ln[Z(T)], nor in fact the partition func- 
tion itself is accessible directly through importance sam- 
pling. Previous calculations of Sa were obtained using 
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FIG. 1: Boundary conditions of the modified simulation, re- 
lated to Equation [6] States with the same labels represent 
where there periodic boundary conditions of the simulation 
cell are, while lines represent the paths of up spins (or bosons, 
depending on context) in the imaginary time direction. 



a thermodynamic integration of internal energies ob- 
tained over many separate simulations (or one "anneal- 
ing" simulation) [T5J [20] or by simulations that sample the 
ratio of two distinct partition functions at a fixed tem- 
perature [21j . However, techniques involving extended 
ensembles have been employed for decades in the calcu- 
lation of free energies and related quantities in a variety 
of Monte Carlo techniques. Below, we discuss the adap- 
tation of one of these, Wang-Landau sampling [22 , to 
the calculation of Renyi entropies in finite-temperature 
QMC techniques. 



B. Wang-Landu sampling 

In contrast to importance sampling where we sam- 
ple states using the thermal weights of the system and 
take an unweighted average of observables, Wang-Landau 
sampling seeks to use the sampling process itself to de- 
termine the distribution of an unknown function. If we 
consider Eq. (|3]) the partition function can be re-written 
as 



Z=^/3"g(n), 



(7) 



where n is the number of operators in the SSE formalism, 
Eq. (§. Since in the SSE the energy can be calculated 
using the estimator [34j E = —(ri^/P, by knowing g{n) 
we are also able to calculate the energy and related vari- 



ables like the specific heat for all temperatures (assuming 
complete knowledge g{n)). In addition, we are able to 
explicitly calculate the partition function for any temper- 
ature (up to some cutoff), something not possible when 
using importance sampling. 

Unlike importance sampling in which we sample states 
with weights derived from the Boltzmann distribution, 
the Wang-Landau method samples using the weighting 
W — g^^{E), (or g"^(n) in the quantum case) ensuring 
that all energies are sampled with equal frequency. Since 
we do not know g{E) a priori we start with an estimated 
density of states, typically flat, and sample the system 
assuming that our estimate is good. If g{E) is correct for 
all energy (up to a scaling factor) then the histogram of 
energies sampled will be flat — if the histogram is larger 
for some energy E' , than our estimate for g{E') is too low 
for that energy, and vice versa. In this way we use the 
sampling histogram to iteratively update the density of 
states until the histogram is flat and g{E) is converged. 



C. Quantum Monte Carlo and Wang-Landau 

In this section we will describe all the changes added to 
the SSE QMC necessary for performing quantum Wang- 
Landau. As mentioned before, previous work |26j has 
successfully implemented Wang-Landau sampling in SSE 
QMC, as well as extensions of the method [55]. Our 
implementation mainly focuses on using Wang-Landau 
(as a proof of concept) and includes recent advances made 
in the classical Wang-Landau formalism [?7] . 

The first challenge with all Wang-Landau techniques 
is representing the sampling function, here g{n), in a 
tractable way. This problem of internal representation 
is that g{n) ranges over 2000 orders of magnitude (for 
larger systems of typical Hamiltonians). Given computer 
precision, the common solution that works here is to only 
keep track of the natural logarithm of g{n). This repre- 
sentation is entirely compatible with the simulation, and 
allows for the largest amount of precision without using 
special data types. For the purpose of the below discus- 
sion, we will let G{n) = \n{g{n)). 

The biggest change from the importance sampling sim- 
ulation is a modification of the diagonal update .26; in 
which the number of operators (the n in g{n)) is changed. 
In the normal simulation, we compare the the weight be- 
fore and after, giving us the probability of adding an 
operator as 
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where i?i,f, is a diagonal operator acting on bond h and 
n is the number of operators. When using Wang-Landau 
this update changes to 



P =1 
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(9) 
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using the current internal estimate of G{n). Then, 
whether we accept or reject the move, we update 



din) 
G{n) 



=gold(f^) X /, 
=GoM{n) + F, 



(10) 
(11) 



using the current n after the update. F = log(/) is 
the refinement parameter used to converge the function 
G{n). The refinement parameter starts as some large 
value (in our case e, the natural log) and when the current 
histogram is found to be fiat (±10 percent of the mean 
for each bin) F is reduced by 



F = Fold/2, 



(12) 



and the histogram is reset. The directed loop update is 
carried out as before [23] (after each diagonal update), 
but G{n) is not changed or used in this update since the 
number of operators does not change. 

The above algorithm works well when the estimate 
of G{n) is far from convergence, but more recent work 
on the classical Wang-Landau algorithm has shown 
that a small saturated error may not converge with this 
method. The concept is simple: when the system is close 
to converged, say the maximum deviation of G{n) from 
the true function is less than our flatness tolerance, then 
the simulation will assume the histogram is flat (within 
tolerance) after a small amount of sampling. This re- 
sults in the reflnement factor getting small very quickly 
without changing G(n) significantly, causing early noise 
in the sampling of G{n) to become frozen in. To correct 
for this, we keep track of the number of full diagonal up- 
dates (update over all elements of the operator list) and 
directed loop updates as a number t. We use the afore- 
mentioned update until we find that F < 1/t, after which 
we ignore the histograms and instead set the refinement 
parameter to 



F = l/t. 



(13) 



This slower refinement is proven to converge the function 
G(n). This choice of F ensures both that F is decreasing, 
and that any saturated error can be overcome through 



enough simulation, as ^^"t ^/n > C can be satisfied 
for any C (representing the saturated error) and t (rep- 
resenting the current Monte Carlo step) given a large 
enough Nmc (representing the number of total Monte 
Carlo steps). 



D. Ensemble optimization 

If we move to more general broad histogram tech- 
niques, a good extension for SSE QMC is the "mini- 
mum round trip" approach |28j . This technique tends 
to reduce the statistical error (shown in Figure [2]) by 
re-shifting the sampling effort to the areas contributing 
most to the error, thereby being more efficient overall. 

To briefly summarize the technique, it keeps track of 
separate histograms for upward and downward moving 



walkers, where the walkers are moving in the one di- 
mensional space of number of operators, n. A walker 
is "upward" if it most recently visited n — 0, and "down- 
ward" walkers arc those that most recently visited n — N 
(where N is our cutoff in the number of operators) . With 
these two histograms we have both the total sampling 
(upward plus downward) plus the rate of change in the 
fraction of upward (or downward) walkers as a function 
of n. The concept is, if there is a point in the simula- 
tion where the fraction of upward walkers change sharply, 
this is interpreted as the simulation having trouble tun- 
nelling through a barrier in configuration space. To pass 
through such barriers more easily, we weight the simu- 
lation to locations where the change in the fraction of 
walkers is large. In addition, the sampling is weighted by 
the inverse number of total samples — this alone would 
simply be Wang-Landau sampling (in that the weights 
would converge to the inverse density of states), but the 
addition of the second condition skews the weights. This 
process is iterated until the modified weights converge. 

With the modified distribution we can sample and take 
the product of the histogram of visits and our skewed 
weights to get the density of states (in Wang-Landau the 
histogram of visits is flat, and so this product is not neces- 
sary in principle), and once we have this we can proceed 
with an identical analysis to our treatment of generat- 
ing the density of states with Wang-Landau alone. We 
implement this method for the 2D transverse field Ising 



model and present and discuss the results in Section III 



E. Calculating observables 

Unlike importance sampling where the calculation of 
observables is aggregated over the course of the simula- 
tion, Wang-Landau sampling produces a function from 
which observables can be calculated after the simulation. 
In this section we will discuss the challenges associated 
with using the function G(n) to calculate observables, 
and the variety of data that we calculate in relation to 
entanglement. 

The first step in using G(n) is normalizing the function, 
since Wang-Landau converges G{a) — G(a -I- 1) for all a, 
but not the absolute magnitude of either. By examining 
the partition function at infinite temperature we see that 



Z(/3 = 0)-^/3"5H=5(0), 



?i=0 



ln[Z(/3 = 0)] = 5(/3 = 0) = G(0), 



(14) 
(15) 



where S'(/3 — 0) is the entropy at infinite temperature. 
In a simulation employing spin-1/2 particles, at infinite 
temperature each spin is independent and the entropy is 



5(/3 = 0) = iV,ln(2), 



(16) 



with N g as the number of spins. In the modified simula- 
tion with n layers, spins not in the region are duplicated 
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(n— 1) times, and act as independent degrees of freedom. 
In this (more general) case, the entropy is 

^„(A,/3 = 0) = (A^, + (a-l)(iV, -iV^))ln(2), (17) 

where Na is the number of spins in region A. Using these 
formulae we can normalize G{n) for any of the simula- 
tions by subtracting a constant amount from all elements 
of G{n) such that G(0) has the correct value. 

In addition, if we know the energy at infinite temper- 
ature we can gather information about G(l). This is 
because if we look at the calculation of energy it can be 
written as 



E : 



lim E 



„G(n)^n-l 



,G(1)-G(0) 



+ c, 



(18) 
(19) 



where C is the sum of all constants subtracted from the 
Hamiltonian in the SSE formalism. This gives us a sec- 
ond check for any model in which the high temperature 
limit of the energy is known. 

Using Eq. ([t]) and the normalized G{n) we get the 
equation for the partition function. 



ln[Z(/3)]=ln[^/3«f 



,G(n) 



(20) 



Although such a calculation does not pose any problems 
analytically, computing such a large term using fixed pre- 
cision data types is a small challenge. To calculate it, first 
we use a second function G'{f3, n) defined as 



G'{(3,n) = G{n) + nln{l3). 



(21) 



Using this function the log of the partition function can 
be written as 



ln[Z(/3)]=ln[^, 



.G'{f!,n) 



(22) 



Finally, we factor the entire equation by K 
max{G"(/3, n)} to get 



HZ{I3)] = ln[K] +ln[^( 



,G'{l3,n)-K 



(23) 



By doing this, the series can be re-exponentiated, 
summed, and the natural logarithm taken of the result 
without the risk of numerical overflow. Other calcula- 
tions, such as the energy, use a similar trick to prevent 
the very large values of the partition function from caus- 
ing problems in computation. 

The method so far allows calculation of estimators de- 
rived from G(n), but gives no estimate on the quality of 
the data. To determine the method's accuracy (without 
comparing to known results) we use the fact that even 
using the 1/t optimization, the initial noise of the sim- 
ulation is only proven to asymptotically converge. By 
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FIG. 2: Comparison of the first derivative of din) using 
Go(n) (arbitrary) as a baseline for i = 0..15. Each color 
represents a different independent simulation of the 2D XXZ 
model of a 32 X 32 lattice. The noise is evenly distributed 
with a larger width for the minimal and maximal values of n. 



using multiple independent simulations we can estimate 
any observable F[G{n)] as 



{F[G[n) 
A(F[G(n) 



[-5:(F[G.(n)]) 

i=l 



(24) 
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(25) 



where (• • • ) represents the statistical average of a quan- 
tity over the simulation and Gi{n) represents the func- 
tion generated from the i"^ simulation. The standard 
deviation of the mean of {F[G{n)]) is represented by 
A(F[G(n)]). We denote the estimator as {F[G{n)]) to 
distinguish our statistical estimate for the observable to 
the analytically exact F[G{n)] if we had precise knowl- 
edge of G(n). 

The noise in G(n) might make it tempting to estimate 
F[G{n)] as 



{F[G{n)]) =F[(G(n))], 
1 ^ 

(Gin)) 



(26) 
(27) 



We have tested that calculating the estimate of the 
observables by Eq. ( 24 ) gives very similar results to 
Eq. ( 26 ) , but it does not lend itself easily to calculating 
the confidence of the the estimate. Comparison of G(n) 
from different simulations show that the deviation from 
the mean is highly correlated — that is to say if Gi{n) is 
larger than average that Gi{n + k) is also likely to be 
larger than average for small k. This error is clarified 
by examining the distribution discrete first derivative of 
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Gi{n) in Fig. [2] That the first derivative is the quan- 
tity that has a random uncorrelated distribution rather 
than the function itself can be understood by under- 
standing that it is the ratio between adjacent elements, 
g{n + l)/ g{n) or G{n + 1) — G{n) , that the Wang-Landau 
simulation converges for all n, and through this the algo- 
rithm is able to reconstruct G{n). 



III. RESULTS FOR MUTUAL INFORMATION 

The Wang-Landau technique outlined above is particu- 
larly suited for studying finite-temperature properties of 
the Renyi entanglement entropies. It's important to note 
that, at finite temperature, the Renyi entropies in general 
no longer obey the property that Sa(pA) — Saips), since 
each quantity picks up volume contributions to its scal- 
ing from thermal mixing as the temperature is increased 
from T = 0. The analogous quantity that is studied at 
finite temperatures is therefore the mutual information 
(MI), 



SaiPA) + SaiPs) — SaipAUs), 



(28) 



which is designed to cancel the volume contributions af- 
fecting each Renyi entropy at T > 0, resulting in a sym- 
metric quantity that reduces to la — '^SaipA) — '^Saips) 
at T = assuming Sa (p) = at zero temperature (which 
is not true if there is ground state degeneracy). In QMC, 
the MI can be calculated for a given region A and its 
compliment (B) using three simulations — two if A and 
B are congruent. In previous work |15j the calculation 
for the Renyi entropy was done using thermodynamic in- 
tegration, extending Eq. (17) to finite temperature by 



integrating energy, similar to how one would calculate 
thermodynamic entropy in classical systems: 

SaiPA) ^-Sc^iA, 13 = 0)+aSil3 = 0) 

+ / {E)j^p,dp'^a / {E)p,dp', (29) 
Jo ' Jo 

where {E) p, {{E) p,) is the energy of the modified simu- 
lation cell (normal simulation) at inverse temperature /?'. 
Using this integration technique, the MI was studied us- 
ing QMC in the context of the anisotropic Heisenberg (or 
XXZ) model HHj- To calculate results using Eq. ([29| sim- 
ulations were run at many temperatures to provide the 
value of {E) ^, over a sufficiently wide range and fine 
temperature mesh to perform numerical integration. By 
studying the finite-size scaling of the a-Renyi entropies, 
the authors discovered that any critical behavior manifest 
at a finite-temperature phase transition Tc is manifest 
as approximate crossings of the MI (for different lattice 
sizes) at and aT^ gD]. 

We first examine the behavior of the MI in the 2D 
square-lattice XXZ model. 



(30) 
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FIG. 3: I2/I, where / is the boundary length, for the 2D 
Heisenberg model with L = 12, comparing the integration 
of importance sampling and Wang-Landau. Notice the sharp 
anomalous behavior in the Wang-Landau for temperatures 
slightly below the cutoff, here chosen to be T = 0.2. This 
means we can only guarantee results down to that tempera- 
ture, but not that they necessarily diverge there. 



using our Wang-Landau algorithm of the previous sec- 
tion, and compare its performance to the conventional 
integration technique [20]. In addition to the 2D 
anisotropic Heisenberg model of Ref. |2D] (A = 4), we 
examine the 2D Heisenberg model(A = 1). To test an ex- 
tended Wang-Landau method, we use the 2D transverse 
field Ising model |35j in the regime where the transverse 
field is weak and it has a finite temperature phase to an 
Ising-like ground state. The Hamiltonian for this model 
is typically written as 



H 



(31) 



where and are Pauli spin operators with eigenval- 
ues ±1, and we study the model at F = J = 1. Finally, 
we examine the 3D Heisenberg model on a cubic lattice 
[36], which also has a finite-temperature transition. 

To begin, we examine the dependence of the MI on 
the cutoff error due to our knowledge of g(n) only up to 
some large (finite) value of n. In the usual SSE QMC 
the number of operators n sampled is fluctuating, but 
for a given j3 we can associate an upper bound on the 
number of operators needed to faithfully represent the 
system, N. This N can be used to choose an appropri- 
ate cutoff in the Wang-Landau algorithm. Choosing to 
have a hard cutoff in n corresponds to a cutoff in reliable 
data at some /?. This can be chosen to be slightly larger 
than the N identified for the cutoff /3 in order to ensure 
that data is still accurate at the f3 of interest — typically 
1.2 to 1.4 times N. Results generated for temperatures 
above the cutoff are accurate, but those below the cutoff 
exhibit anomalous behavior, as effectively g{n) = for 
n > N (or equivalently, G{n) — —00). The result of at- 
tempting to generate results below this cutoff is shown 
in Fig. |3] G{n) is a smooth function for large enough n. 



7 




4 6 

Temperature 



FIG. 4: The energy per spin and mutual information divided 
by boundary length for the 2D XXZ model with A = 4. No- 
tice the crossings of the mutual information very close to Tc 
and 2Tc, marked with vertical lines. Different lines correspond 
to lattices are of size 4, 6, 8, 10, 12, 14, 16, 18, 20 and 22 (from 
lowest to highest at peak). 
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FIG. 5: The energy per spin and mutual information divided 
by system size for the 2D Heisenberg model. This model 
shows no fixed crossing in 72 /Z. Neighboring sizes cross at 
progressively lower temperatures. Different lines correspond 
to lattices are of size 4, 6, 8, 10, 12, 14, 16, 18, 20 and 22 (from 
lowest to highest central peak). 



SO there is room for methods which generate data using 
analytic extensions of G{n) to reduce the anomalous be- 
havior near the lower temperature bound, but in general 
exactly knowing such an analytic form (for all n) would 
be equivalent to having an analytical expression for the 
energy as a function of temperature. 

We now briefly examine the finite-size scaling of the 
MI for several examples of A in the XXZ model. To be- 
gin, recall that the anisotropic Heisenberg model (with 
A = 4) undergoes a phase transition to a z-axis polarized 
antiferromagnetic ground state as the temperature is low- 
ered. The transition temperature was determined previ- 
ously using the fourth order binder cumulant of the stag- 
gered magnetization, finding Tc — 2.42 [20] ■ The energy 
per spin and mutual information I2 divided by boundary 
size I = 2L is shown in Fig. |4] The mutual information 
per system size for different L exhibits a universal scal- 
ing |20] manifest as approximate crossings at Tc, and 2Tc 
(when using I2 from the second Renyi entropy). The in- 
set in Figure |4] shows a closeup of the crossing, with the 
uncertainty of the data points calculated from Eq. (25). 



It is well known that the 2D isotropic Heisenberg model 
(A = 1) has no finite temperature phase transition. 
When examining the mutual information, we do not see 
any fixed crossing in the curves as a function of L, as 
shown in Fig.jSj rather, the crossings move towards T = 
for increasing L. 

To test the extended broad histogram sampling algo- 
rithm |28j we illustrate results from another Hamiltonian, 
the transverse field Ising model. On the 2D square lattice 




3 4 5 6 
Temperature 



FIG. 6: The mutual information divided by system size for 
the 2D transverse field Ising model using the extended broad 
sampling approach at F = J = 1. Inset shows a closeup of 
the lower crossing of mutual information at around T = 2.21, 
crossing which we expect to occur at the transition tempera- 
ture. Different lines correspond to lattices are of size 6, 8, 10 
and 12 (from lowest to highest central peak). 



the classical Ising model has a transition temperature of 
Tc = 2J/ln(l + 72) w 2. 269 J from Onsager's exact re- 
sult [37] . When we add a small transverse field we expect 
that the transition temperature will be slightly lowered, 
and as the critical point of this model for the two di- 
mensional system is Fc = 3.044J we are well within the 
Ising like phase at F = J [38] . Figure |6] shows the mu- 
tual information curves for this model for a small set of 
sizes, showing crossings around Tc ~ 2.21, slightly below 
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FIG. 7: The energy per spirr and mutual information divided 
by surface area (still denoted /) for the 3D Heisenberg model. 
The crossings of the mutual information very close to Tc and 
2Tc, marked with vertical lines. Different lines correspond to 
lattices are of size 4, 6, 8 and 10 (from lowest to highest central 
peak). 

the non-perturbed result, as expected. The inset shows 
a very tight error bound for the data points when using 
the extended broad histogram technique. 

Finally, we examine the performance of this method 
in higher dimensions using the 3D isotropic Heisenberg 
model. Unlike the 2D case, the 3D Heisenberg model 
does have a finite temperature phase transition, occur- 
ring for A = 1 at Tc = 0.946 [36J. In this case, the 
regions A and B used to calculate the Renyi entangle- 
ment entropy consisted of a periodic plane spanning the 
3D simulation cell, cutting it into two halves. MI results 
are illustrated in Fig. [7] for a largest size of 10 x 10 x 10. 
In Fig. [7j the MI per surface area of this plane exhibits 
crossings approach the critical temperature within error 
bars. Interestingly, finite size effects are magnified near 
2Tc, when compared to the 2D transition studied above. 
It would be interesting to develop a full scaling theory for 
the 3D Heisenberg transition to examine these effects in 
the MI, similar to what was done previously in 2D [5D]. 

IV. DISCUSSION 

Although Wang-Landau simulations in themselves are 
not new [22] , the use of such non-canonical sampling tech- 
niques in quantum Monte Carlo algorithms is still being 
explored. In general, techniques like Wang-Landau tend 
to suffer from slower convergence for a specific parameter 
(i.e. the energy at a specific temperature) when compared 
to importance sampling techniques, but the trade-off is 
that Wang-Landau generates this result for all tempera- 
tures, something that would require many more simula- 
tions using importance sampling. In the end, the decision 
of the type of sampling to use should depend on the na- 
ture of the most important observable. 



Wang-Landau can also be extended to measure any 
quantity that could be measured using traditional SSE. 
Instead of just simply having g{n), we now associate any 
number of other observables (like magnetization) with 
each g{n). By looking at the average value sampled for 
a given n, when we re-weight all the g{n) we can cal- 
culate these quantities for any temperature. One dif- 
ference between Wang-Landau and importance sampling 
here is that because all energies are sampled, the simula- 
tion does not tend to get stuck in local configurations. A 
system with a symmetry broken ground state, like the 2D 
Ising model, will sample all ground states with higher fre- 
quency than most importance sampling methods. This 
requires some modified sampling styles similar to that 
needed with parallel tempering, but it is not difficult to 
implement. 

We have demonstrated that the Wang-Landau tech- 
nique can be successfully used along with a modified SSE 
formalism [15] to generate the Renyi entanglement en- 
tropies, and associated mutual information, over a wide 
range of temperature. We confirm that extending simu- 
lations to use the modified 1/t schedule ensures that any 
simulation error is systematically reduced by continued 
sampling. Improvements using optimized ensembles al- 
lows for a similar analysis once the density of states is 
extracted. 

Such broad histogram techniques can also be extended 
to work in other parameters of the Hamiltonian, such as 
the size of the entangled region A in the replica trick. A 
Wang-Landau in this parameter would allow the calcula- 
tion of the entanglement entropy at a fixed temperature 
without integration. In fact, the implementation by Hu- 
meniuk |21j is identical to the implementation of such a 
Wang-Landau in region size, with only minor changes. 
As the accumulation of error is still a concern, it is not 
trivially clear that integrating in energy or integrating 
in region size leads to a smaller error for all cases. This 
should be explored in the future. 

Using our Wang-Landau method, we have obtained 
mutual information results for the 2D anisotropic Heisen- 
berg model, reproducing the crossings at the critical tem- 
perature [5D]. For the isotropic case, no critical scaling of 
the mutual information is observed, consistent with the 
lack of a finite-temperature phase transition in the model. 
We have also implemented a QMC for the 2D transverse 
field Ising model, where the mutual information again 
shows evidence of the expected finite temperature phase 
transition. Finally, we examine the mutual information 
in a three dimensional Heisenberg model, and show that 
crossings occur at Tc and approximately 2Tc. 
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